%% The actual Travel model and modifications
%
% Input:  PreferenceFactor (used for minimization),cfg,exvars,depvars,indepvars,r,t
% Output: Residuals (depeding on calibration type; only during calibrations)
%         depvars structure
%
function varargout = ModelTravelOptimizedTW(cfg,exvars,depvars,indepvars,r,t,varargin)
%% Calls the function in a method that allows the model to minimize the difference betweeen TTB and actual use
if isempty(varargin)
   [depvars,indepvars] = SelectSolver(cfg,exvars,depvars,indepvars,r,t,['''depvars.TimeWeight''',',x(1)'],1);
   [varargout{1}, varargout{2}] = RunModel(cfg,exvars,depvars,indepvars,r,t,varargin{:});
else
   varargout{1} = RunModel(cfg,exvars,depvars,indepvars,r,t,varargin{:});
end

end

function varargout = RunModel(cfg,exvars,depvars,indepvars,r,t,varargin)

    for n = 2:2:length(varargin)
         eval( [varargin{n-1} '(r,time(1),:)=varargin{n};']);   
    end
    
    modes                               = indepvars.UseModes(r,1,:); % get rid of modes not used by Schafer (e.g. High speed train)    
    depvars.PreferenceFactor(r,:,:)   = repmat(indepvars.PreferenceFactorStart(r,:,:), [1 length(cfg.Time) 1]) ...
            + repmat(cfg.Time,size(indepvars.PreferenceFactorStart(r,:,:))) .* repmat(indepvars.PreferenceFactorSlope(r,:,:),[1 length(cfg.Time) 1]);

for t = time
%% Girod et. al TIMER model
    if sum(strcmp('depvars.TimeWeight',varargin))==0
    depvars.TimeWeight(r,t) =...
        depvars.TimeWeight(r,t-1) .* ...
        sum( depvars.TravelDemand(r,t-1,modes )./exvars.Speed(r,t-1,modes),3) ...
        ./ indepvars.TTB(r);
    end

    depvars.RelativeModeCostTravel(r,t,modes ) =...
    exvars.TravelPrice(r,t,modes ) ./ (mean(exvars.TravelPrice(r,t,modes )) .* depvars.PreferenceFactor(r,t,modes))...
    + 1./exvars.Speed(r,t,modes) ./ mean(1./exvars.Speed(r,t,modes)) .* depvars.TimeWeight(r,t);

    depvars.ModeSplitOptimum(r,t,modes ) =...
    exp(-indepvars.Lambda(r) .* depvars.RelativeModeCostTravel(r,t,modes )) ...
    / sum(exp(-indepvars.Lambda(r) .* depvars.RelativeModeCostTravel(r,t,modes )));
    
    depvars.TMB(r,t)             = (0.035 + 0.085 * sum(depvars.ModeSplit(r,t,5:end),3)) * indepvars.TMBMultiplier(r); 

    depvars.TravelDemand(r,t,modes) = ...
       depvars.ModeSplit(r,t,modes ) .* depvars.TMB(r,t) .* exvars.Income (r,t) ./ sum(depvars.ModeSplit(r,t,modes ) .* exvars.TravelPrice(r,t,modes),3);
    
end

%% Write output
    varargout{1} = mean(depvars.ObjectiveFunction(r,time));
if nargout==2
    varargout{2} = depvars;
end
    
%% Visualizes minimization process
if cfg.LiveMinimizeVisuals{1} > 1 && ismember(time(1),cfg.LiveMinimizeVisuals{2}) && length(depvars.Resid)<cfg.LiveMinimizeVisuals{1}; % bit clumsy to make the code faster
     GeneralPlotting('LiveMinimizeVisuals','Dont save',cfg,exvars,depvars,indepvars,r,time);
elseif cfg.PlotMinimize == 1
    GeneralPlotting('SaveToGlobal',' ',cfg,exvars,depvars,indepvars,r,time);
end
%% Determine the errors
   
    depvars.ModeWeight(r,t,modes) = sqrt(exvars.TravelDemand(r,t,modes))./sum(sqrt(exvars.TravelDemand(r,t,modes)),3);% Ensures not all modes are equal, but that car is not the only important mode
    
    if ismember('depvars.TimeWeight',varargin(1:2:end)) && ~ismember('depvars.PreferenceFactor',varargin(1:2:end))
    depvars.ObjectiveFunction(r,t) = (sum(depvars.TravelDemand(r,t,modes)./exvars.Speed(r,t,modes)) - ...
    indepvars.TTB(r,t)).^2./indepvars.TTB(r,t)^2;           
    else
    if t<36 % ensures function is only called with historic data
    depvars.ObjectiveFunction(r,t) = sum((depvars.TravelDemand(r,t,modes) - exvars.TravelDemand(r,t,modes)).^2 ...
    ./ exvars.TravelDemand(r,t,modes).^2 .* depvars.ModeWeight(r,t,modes).^2); % used for optimization
    depvars.Resid(r,t)            = sqrt(sum((sum(depvars.TravelDemand(r,t,modes),3) - sum(exvars.TravelDemand(r,t,modes),3) ...
    ./ sum(exvars.TravelDemand(r,t,modes),3)).^2,2)./length(time)); % used to denote the error (NRMSE)
    end
    end

%% Write output
    varargout{1} = mean(depvars.ObjectiveFunction(r,t));
if nargout==2
    varargout{2} = depvars;
end
    
%% Visualizes minimization process
if cfg.LiveMinimizeVisuals{1} > 1 && ismember(time(1),cfg.LiveMinimizeVisuals{2}) && length(depvars.Resid)<cfg.LiveMinimizeVisuals{1}; % bit clumsy to make the code faster
     GeneralPlotting('LiveMinimizeVisuals','Dont save',cfg,exvars,depvars,indepvars,r,time);
elseif cfg.PlotMinimize == 1
    GeneralPlotting('SaveToGlobal',' ',cfg,exvars,depvars,indepvars,r,time);
end
    
end